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Abstract 

This  paper  addresses  the  development  of  a  unified  framework  for  quantifying  hysteresis  and  con¬ 
stitutive  nonlinearities  inherent  to  ferroelectric,  ferromagnetic  and  ferroelastic  materials.  Because 
the  mechanisms  which  produce  hysteresis  vary  substantially  at  the  microscopic  level,  it  is  more  nat¬ 
ural  to  initiate  model  development  at  the  mesoscopic,  or  lattice,  level  where  the  materials  share 
common  energy  properties  along  with  analogous  domain  structures.  In  the  first  step  of  the  model 
development,  Helmholtz  and  Gibbs  energy  relations  are  combined  with  Boltzmann  theory  to  con¬ 
struct  mesoscopic  models  which  quantify  the  local  average  polarization,  magnetization  and  strains 
in  ferroelectric,  ferromagnetic  and  ferroelastic  materials.  In  the  second  step  of  the  development, 
stochastic  homogenization  techniques  are  invoked  to  construct  unified  macroscopic  models  for  non- 
homogeneous,  polycrystalline  compounds  exhibiting  nonunifornr  effective  fields.  The  combination  of 
energy  analysis  and  homogenization  techniques  produces  low-order  models  in  which  a  number  of  pa¬ 
rameters  can  be  correlated  with  physical  attributes  of  measured  data.  Furthermore,  the  development 
of  a  unified  modeling  framework  applicable  to  a  broad  range  of  ferroic  compounds  facilitates  ma¬ 
terial  characterization,  transducer  development,  and  model-based  control  design.  Attributes  of  the 
models  are  illustrated  through  comparison  with  piezoceramic,  magnetostrictive  and  shape  memory 
alloy  data  and  prediction  of  material  behavior. 
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1  Introduction 


Ferroelectric,  ferromagnetic  and  ferroelastic  materials  share  a  number  of  mesoscopic  and  macroscopic 
attributes  including  the  formation  of  analogous  domain  structures  and  the  presence  of  hysteresis  and 
constitutive  nonlinearities  as  illustrated  in  Figure  1.  While  hysteresis  and  constitutive  nonlinearities 
are  inherent  to  all  presently  employed  ferroic  compounds,  the  degree  and  severity  of  these  effects 
can  often  be  mitigated  by  restricting  drive  levels,  employing  appropriate  drive  electronics,  or  in¬ 
corporating  feedback  loops  in  transducers.  In  many  applications,  however,  hysteresis  and  material 
nonlinearities  are  unavoidable  and  must  be  incorporated  in  models  and  subsequent  control  designs  to 
achieve  the  full  capabilities  of  the  materials.  To  illustrate,  we  consider  the  behavior  of  certain  repre¬ 
sentative  ferroelectric,  ferromagnetic  and  ferroelastic  materials  employed  in  present  high  performance 
transducer  applications. 

Piezoceramic  actuators  are  employed  for  nanoscale  positioning  due  to  their  high  set  point  accu¬ 
racy  and  broadband  capability.  In  present  atomic  force  microscope  (AFM)  or  scanning  tunneling 
microscopic  (STM)  designs,  PID  or  robust  control  laws  are  employed  to  mitigate  hysteresis  in  the 
relation  between  input  fields  or  voltages  and  strains  generated  by  the  device  [7,  28,  87,  88,  89,  90]. 


Strain 

(c) 

Figure  1.  Hysteresis  and  constitutive  nonlinearities  exhibited  by  various  ferroic  compounds, 
(a)  PZT5A  data,  (b)  Terfenol-D  data,  and  (c)  NiTi  data. 


1 


This  is  effective  for  low  scan  rates  and  has  led  to  the  phenomenal  success  of  the  devices.  However,  at 
the  higher  scan  rates  required  for  mass  product  diagnostics  (e.g.,  quality  evaluation  of  semiconductor 
chips)  or  real-time  monitoring  of  biological  processes  (e.g.,  observation  of  protein  folding  dynamics), 
noise  to  data  ratios  and  diminishing  high  pass  characteristics  of  control  filters  preclude  the  sole  use 
of  feedback  control  laws  to  mitigate  hysteresis,  and  control  designs  utilizing  model-based  inverse 
compensators  are  under  investigation  [24,  25,  111,  124].  In  other  applications  utilizing  piezoce¬ 
ramic  actuators,  the  use  of  charge  or  current  controlled  amplifiers  can  essentially  eliminate  hysteresis 
[64,  65,  66,  67,  68].  However,  this  mode  of  operation  can  be  prohibitively  expensive  as  compared 
with  the  more  commonly  employed  voltage  control  amplifiers,  and  current  control  is  ineffective  if 
maintaining  DC  offsets  as  required  for  numerous  applications  -  e.g.,  the  x-stage  in  an  AFM  must 
hold  a  specified  position  while  a  sweep  is  performed  in  the  y-stage. 

In  magnetic  materials,  significant  research  is  focused  on  the  development  of  compositions  and 
compounds  which  maximize  performance  specifications  such  as  output  strains  and  blocked  forces 
while  minimizing  hysteresis.  This  has  motivated  the  development  of  new  rare  earth  transducer  com¬ 
pounds,  such  as  Galfenol  [21,  22],  as  well  as  ferromagnetic  shape  memory  alloy  (FSMA)  compounds 
[20,  37,  54,  78,  116].  However,  the  materials  with  the  highest  work  density  metrics  still  exhibit  signif¬ 
icant  hysteresis  and  constitutive  nonlinearities  which  must  be  accommodated  in  transducer  design. 
Moreover,  for  other  technologies  including  magnetic  recording,  a  high  degree  of  hysteresis  is  actually 
required  to  provide  the  history  required  for  the  application. 

In  the  class  of  ferroelastic  materials,  shape  memory  alloys  (SMA)  are  being  increasingly  considered 
for  civil,  aeronautic,  aerospace  and  industrial  applications  which  require  significant  passive  damping 
or  utilize  the  high  work  output  densities  exhibited  by  the  materials.  Because  the  energy  dissipated 
by  ferroic  materials  is  proportional  to  the  area  of  the  hysteresis  loop,  pseudoelastic  operating  regimes 
which  maximize  hysteresis  are  required  when  employing  SMA  as  tendons  to  attenuate  earthquake 
or  wind-induced  vibrations  in  buildings  or  as  fibers  to  eliminate  vibrations  in  articulated  antennas 
or  membrane  mirrors  -  see  [4,  33,  34,  35,  121]  for  recent  civil  applications  and  [91,  92]  for  details 
regarding  the  modeling  of  hysteresis-induced  damping  behavior.  The  utilization  of  temperature- 
induced  phase  transitions  to  provide  actuator  capabilities  is  under  intense  investigation  in  the  context 
of  microelectromechanical  systems  (MEMS),  thin  film  SMAs,  and  nricroactuator  applications  since 
surface  area  to  volume  ratios  in  these  geometries  promote  rapid  cooling  and  hence  high  frequency 
drive  capabilities.  For  example,  Ho  et  al.  report  on  a  thin  film  micropump  capable  of  operating  at 
300  Hz  [46].  Additional  details  regarding  thin  film  SMA  applications  can  be  found  in  [39,  61,  69]  and 
discussion  regarding  the  development  of  SMA-based  microgrippers  relevant  for  microrobotics  and 
microassembly  is  provided  in  [59,  84,  85].  To  achieve  both  the  stress-induced  damping  capabilities 
and  temperature-induced  actuator  capabilities  of  SMA  devices,  it  is  necessary  to  quantify  the  inherent 
hysteresis  and  constitutive  nonlinearities  in  a  manner  feasible  for  transducer  design  and  real-time 
control  implementation. 

Finally,  there  is  an  increased  focus  on  the  design  of  hybrid  transducers  which  utilize  comple¬ 
mentary  properties  of  constituent  materials.  For  example,  the  drive  characteristics  and  90°  phase 
shift  between  Terfenol-D  and  PZT  or  PMN  is  being  utilized  to  construct  hybrid  transducers  hav¬ 
ing  improved  energy  efficiencies  and  frequency  bandwidths  [15,  16,  36,  55].  However,  the  hysteresis 
and  nonlinearities  inherent  to  the  constituent  materials  must  be  accommodated  to  achieve  optimal 
performance  with  hybrid  designs. 

Hence  from  both  the  fundamental  perspective  of  material  characterization  and  future  material  de¬ 
velopment,  and  the  practical  perspective  of  transducer  design  and  model-based  control  development, 
it  is  advantageous  to  develop  unified  frameworks  for  quantifying  hysteresis  and  material  nonlinear¬ 
ities  inherent  to  ferroelectric,  ferromagnetic  and  ferroelastic  materials.  Material  development  and 
characterization  necessitate  a  high  degree  of  accuracy  whereas  the  practical  requirements  of  trans- 
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ducer  design  and  real-time  control  implementation  limit  the  complexity  of  models  and  favor  models 
which  can  be  efficiently  constructed  and  updated  to  accommodate  changing  operating  conditions.  In 
this  paper,  we  employ  a  multiscale  development  to  provide  a  unified  framework  for  quantifying  the 
hysteresis  and  constitutive  nonlinearities  inherent  to  a  broad  range  of  ferroic  compounds.  The  goal 
is  to  provide  sufficient  accuracy  for  material  characterization  while  maintaining  a  complexity  level 
which  facilitates  transducer  and  control  design  in  addition  to  real-time  control  implementation. 

As  summarized  in  Table  1,  ferroelectric,  ferromagnetic  and  ferroelastic  materials  share  a  num¬ 
ber  of  common  traits  including  the  presence  of  multiple  domains,  which  are  separated  by  domain 
boundaries  or  walls,  and  field-induced  phase  transitions  in  the  vicinity  of  the  Curie  temperature  Tc. 
A  unified  analysis  of  these  shared  properties  dates  back  at  least  to  the  work  of  Nye  [77]  and,  in  1970, 
Aizu  presented  a  unified  treatment  of  certain  symmetry  properties  of  ferroelectric,  ferromagnetic  and 
ferroelastic  compounds  which  he  collectively  referred  to  as  ferroic  compounds  [5] .  Shared  attributes 
of  ferroic  compounds  have  subsequently  been  exploited  by  a  number  of  researchers  when  categorizing 
properties  of  the  materials  [76,  109,  120]. 

The  physical  mechanisms  which  produce  hysteresis  and  constitutive  nonlinearities  in  ferroelectric, 
ferromagnetic  and  ferroelastic  compounds  differ  significantly  at  the  microscopic  scales.  In  ferroelec¬ 
tric  materials,  hysteresis  is  partially  attributed  to  dipole  switching  and  domain  wall  losses  whereas 
moment  rotation  and  domain  wall  losses  produce  hysteresis  in  ferromagnetic  compounds.  While  a 
number  of  analogies  between  hysteresis  mechanisms  in  ferroelectric  and  ferromagnetic  compounds 
can  be  made  based  on  shared  domain  properties  of  the  materials,  the  actual  molecular  mechanisms 
and  scales  on  which  domain  properties  occur,  differ  substantially.  Furthermore,  hysteresis  in  shape 
memory  alloys  is  due  to  phase  transitions  between  austenite  and  martensite  variants.  Hence  the 
construction  of  unified  energy-based  models  at  microscopic  scales  is  not  considered  feasible. 

At  the  lattice,  or  mesoscopic,  and  macroscopic  scales,  however,  shared  domain  and  energy  prop¬ 
erties  can  be  utilized  to  construct  unified  models  for  a  broad  range  of  ferroic  materials.  In  the  first 
step  of  the  present  development,  Helmholtz  and  Gibbs  energy  relations  are  derived  at  the  lattice  level 
and  employed  in  Boltzmann  theory  to  obtain  evolution  equations  which  quantify  the  local  average 
polarization,  magnetization  and  strains  in  ferroelectric,  ferromagnetic  and  ferroelastic  materials.  For 
homogeneous,  single  crystal  compounds,  these  local  relations  can  be  extended  throughout  the  mate¬ 
rial  to  obtain  macroscopic  constitutive  relations  which  quantify  the  bulk  material  behavior.  In  the 
second  step  of  the  model  development,  stochastic  homogenization  techniques  are  employed  to  incor¬ 
porate  material  nonhomogeneities,  polycrystallinity,  and  nonuniform  effective  fields  found  in  typical 
ferroic  compounds.  The  resulting  macroscopic  models  are  low  order,  and  hence  efficient  to  imple¬ 
ment,  with  parameters  that  can  be  correlated  with  properties  of  the  data.  Furthermore,  the  unified 
models  guarantee  the  closure  of  biased  minor  loops  and  enforce  the  ‘deletion’  property  in  quasistatic 


Ferroelectric 

Ferromagnetic 

Ferroelastic 

Polarization 

Magnetization 

Strain 

Electric  field 

Magnetic  field 

Stress 

Paraelectric  phase 

Paramagnetic  phase 

Austenite  phase 

Ferroelectric  phase 

Ferromagnetic  phase 

Martensite  phase 

Ferroelectric  domain  walls 

Bloch  or  Neel  walls 

Boundaries  between  variants 

Devonshire  theory 

Mean  field  (Weiss)  theory 

Landau  theory 

Micromechanical  ferroelectric  theory 

Micromagnetic  theory 

Ginzburg-Landau  theory 

Table  1.  Analogies  between  physical  properties  of  ferroic  materials. 
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drive  regimes,  and  incorporate  various  temperature  dependencies  and  relaxation  mechanisms.  The 
resulting  unified  modeling  framework,  which  has  its  genesis  in  the  SMA  theory  of  Achenbach  and 
Muller  [1,  2],  thus  quantifies  a  wide  range  of  phenomena  common  to  ferroelectric,  ferromagnetic  and 
ferroelastic  compounds  and  pertinent  to  transducer  designs  which  exploit  the  unique  capabilities  of 
these  materials. 

To  place  this  framework  in  perspective,  it  is  useful  to  compare  it  with  existing  energy-based 
and  phenomenological  methodologies  used  to  quantify  hysteresis  in  ferroic  compounds.  We  refer  the 
reader  to  [18,  48,  51,  52,  62,  86,  108,  110]  for  a  summary  of  nonlinear  and  hysteresis  models  for 
ferroelectric  and  piezoelectric  compounds,  [17,  29,  30,  31,  56,  58,  63,  71]  for  ferromagnetic  materials, 
and  [10,  40,  49,  50,  80,  81,  82,  99]  for  an  overview  of  models  for  shape  memory  alloys.  We  summarize 
next  investigations  focused  on  the  development  of  unified  models  for  ferroic  compounds. 

The  theory  in  [109]  exploits  the  shared  domain  structure  of  ferroic  materials  to  construct  unified 
hysteresis  models.  In  that  development,  statistical  mechanics  tenets  are  invoked  to  first  construct 
anhysteretic  models  for  idealized  materials  devoid  of  pinning  sites.  Electrostatic,  magnetostatic 
and  elastic  energy  relations  are  then  employed  to  quantify  the  reversible  and  irreversible  energy 
required  to  translate  domain  walls  in  the  materials.  The  resulting  models  are  low-order  but  have  the 
disadvantage  that  it  is  difficult  to  enforce  closure  of  biased  minor  loops  without  a  priori  knowledge  of 
turning  points.  We  point  out  that  the  asymptotic  models  derived  in  the  present  framework  through 
statistical  mechanics  principles  are  the  Ising  model  employed  in  the  theory  of  [109]  to  quantify  the 
equilibrium  anhysteretic  behavior  of  ferromagnetic  and  ferroelectric  materials. 

A  variety  of  phenomenological  techniques  have  also  been  developed  to  construct  unified  hys¬ 
teresis  models  for  ferroic  compounds.  In  the  theory  developed  by  Soukhojak  and  Chiang  [115], 
phenomenological  principles  are  used  to  construct  rheological  models  which  quantify  the  time  and 
frequency-dependent  behavior  of  a  range  of  ferroic  materials.  Preisach  techniques  have  also  been 
employed  for  characterizing  hysteresis  and  nonlinearities  in  numerous  ferroic  compounds  due  to  their 
generality  and  rigorous  mathematical  foundation  [3,  8,  11,  41,  71,  60,  83,  119].  The  generality 
of  Preisach  models  also  constitutes  a  weakness  in  the  technique  since  it  is  difficult  to  incorporate 
known  physics  or  properties  of  the  data  when  determining  parameters  in  the  models.  As  discussed 
in  [31,  63],  extensions  to  classical  Preisach  theory  are  also  required  when  incorporating  reversible  ef¬ 
fects  and  temperature  and  frequency  dependencies,  and  relaxing  congruency  properties  to  model  the 
incongruent  behavior  exhibited  by  materials.  In  [112],  it  is  illustrated  that  the  framework  presented 
here  provides  an  energy  basis  for  extended  Preisach  models  with  four  important  differences,  (i)  The 
energy  derivation  yields  physical  parameters  which  can  be  correlated  with  properties  of  the  data  - 
this  facilitates  model  construction  and  updating  to  accommodate  changing  system  inputs,  (ii)  The 
temperature  dependencies  and  relaxation  mechanisms  are  incorporated  in  the  basis,  or  hysterons, 
rather  than  in  the  densities,  measures  or  parameters  -  this  eliminates  the  necessity  of  lookup  tables 
or  vector-valued  parameters  when  employing  the  models  in  variable  temperature  or  certain  dynamic 
regimes,  (iii)  The  model  does  not  erroneously  enforce  congruency  and  hence  does  not  need  to  be 
modified  to  achieve  the  noncongruency  exhibited  in  a  number  of  operating  regimes,  (iv)  The  model 
automatically  incorporates  reversibility  through  the  energy  derivation  of  the  hysterons  and  hence 
does  not  require  modification  to  accommodate  low  drive  regimes. 

In  Section  2,  we  develop  mesoscopic  energy  relations,  differential  equations  quantifying  the  evo¬ 
lution  of  phase  fractions,  and  asymptotic  relations  characterizing  the  average  local  polarization, 
magnetization  and  strains  in  ferroelectric,  ferromagnetic  and  ferroelastic  compounds.  The  theory  is 
extended  in  Section  3  to  nonhomogeneous,  poly  crystalline  compounds  exhibiting  nonunifornr  effec¬ 
tive  fields  by  considering  parameters  such  as  coercive  fields  and  relative  stresses  to  be  manifestations 
of  underlying  distributions.  Attributes  of  the  unified  models  are  then  illustrated  in  Section  4  through 
comparison  with  PZT5A,  Terfenol-D  and  NiTi  data  as  well  as  prediction  of  material  behavior. 
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2  Mesoscopic  Models 


We  summarize  here  the  development  of  mesoscopic  models  derived  through  the  construction  of  ap¬ 
propriate  Helmholtz  and  Gibbs  energy  relations  at  the  lattice  level  with  Boltzmann  theory  used  to 
quantify  the  probability  of  achieving  various  energy  states.  When  combined  with  differential  equa¬ 
tions  which  quantify  temperature  changes  within  the  materials,  this  yields  a  coupled  set  of  evolution 
equations  quantifying  the  local  average  polarization,  magnetization  and  strains  in  ferroelectric,  fer¬ 
romagnetic  and  ferroelastic  compounds.  For  homogeneous,  single  crystal  compounds  having  uniform 
fields  or  stresses,  these  local  relations  apply  throughout  the  material  and  hence  provide  macroscopic 
models  which  characterize  hysteresis  and  constitutive  nonlinearities.  These  local,  mesoscopic  rela¬ 
tions  are  extended  in  Section  3  to  accommodate  the  effects  of  polycrystallinity,  material  defects, 
inclusions,  nonhomogeneities,  and  variable  effective  fields  and  stresses.  To  simplify  the  discussion, 
we  develop  both  the  mesoscopic  and  macroscopic  polarization,  magnetization  and  strain  relations 
in  the  context  of  piezoceramic,  magnetostrictive  and  shape  memory  compounds.  However,  the  the¬ 
ory  is  sufficiently  general  to  encompass  a  broad  range  of  ferroelectric,  ferromagnetic  and  ferroelastic 
materials. 

To  provide  a  general  framework  for  formulating  the  unified  models,  we  consider  an  order  parame¬ 
ter  e  and  external  fields  p  that  are  thermodynamically  conjugate  to  e.  For  ferroelectric,  ferromagnetic 
and  ferroelastic  materials,  natural  choices  for  e  are  respectively  the  polarization  P,  magnetization 
M,  and  strains  e.  The  corresponding  external  fields  p  are  taken  to  be  the  electric  field  E,  magnetic 
field  H,  and  stresses  a.  For  temperatures  T,  we  denote  Helmholtz  energy  relations  by  ip(e,  T)  and 
Gibbs  energy  relations  by  G(e,ip,T). 

In  the  absence  of  applied  fields,  the  Helmholtz  relations  provide  a  natural  measure  of  energy  and, 
under  the  assumption  of  differentiability,  thermodynamic  equilibria  are  provided  by  the  necessary 
condition 
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In  the  presence  of  an  applied  field,  the  Gibbs  relation 
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G(e,ip,T)  =  ^(e,T)  -  pe 
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is  used  to  quantify  the  total  energy,  and  equilibrium  states  are  specified  by  the  necessary  conditions 
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The  relation  (3)  can  be  interpreted  as  providing  conditions  under  which  the  order  parameter  adjusts 
to  balance  the  internal  energy  with  work  due  to  the  applied  field.  As  will  be  noted  in  subsequent 
discussion,  the  conditions  (3)  also  provide  asymptotic  conditions  specifying  the  local  polarization, 
magnetization  and  strains  at  the  lattice  level  as  well  as  macroscopic  constitutive  relations  if  elec¬ 
tromechanical  or  magnetomechanical  coupling  is  included  in  the  Gibbs  energy  relation. 


2.1  Piezoceramic  Materials 

Presently  employed  piezoceramic  compounds  are  comprised  of  PbTii_2,03  and  PbZr.,,0.3  with  x 
chosen  to  optimize  electromechanical  coupling.  To  simplify  the  discussion,  we  focus  on  PbTiC>3 
to  motivate  the  construction  of  appropriate  energy  relations.  At  temperatures  T  above  the  Curie 
temperature  Tc,  PbTiC>3  exhibits  a  cubic  form,  whereas  it  has  the  tetragonal  form  depicted  in 
Figure  2a  for  T  <  Tc.  A  corresponding  potential  or  Helmholtz  energy  profile  associated  with  the 
position  of  the  Ti  cation  along  the  X3  axis  is  depicted  in  Figure  2b.  In  the  presence  of  an  applied 
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Figure  2.  (a)  Tetragonal  structure  of  PbTiC>3  for  T  <TC  and  resulting  spontaneous  polarization, 
(c)  Helmholtz  energy  as  a  function  of  Ti  position  along  the  X3-axis. 


field  E,  the  energy  landscape  distorts  in  the  manner  depicted  in  Figure  3.  In  the  absence  of  thermal 
energy,  the  local  minimum  associated  with  the  stable  equilibrium  disappears  when  the  field  reaches 
the  local  coercive  field  value  Ec  and  it  becomes  energetically  favorable  for  the  Ti  cation  to  change 
configurations.  This  produces  a  dipole  switch  and  subsequent  hysteresis  in  the  relation  between 
E  and  the  local  average  polarization  P.  Whereas  other  mechanisms,  including  thermal  activation, 
contribute  to  the  hysteresis  and  constitutive  nonlinearities  exhibited  by  piezoceranric  materials,  this 
source  of  dipole  switching  helps  motivate  appropriate  formulations  for  the  Helmholtz  and  Gibbs 
energy  relations. 

2.1.1  Helmholtz  and  Gibbs  Energy  Relations 

Following  the  development  in  [114],  we  consider  a  uniform  lattice  of  volume  V  and  mass  v  having 
N  cells.  We  also  make  the  assumption  that  each  cell  has  two  possible  orientations  s%  =  ±1  and 
dipole  moment  p.  We  let  To  denote  the  energy  required  to  reorient  a  single  dipole  when  the  lattice 
is  completely  ordered. 

The  Helmholtz  energy  for  the  lattice  is 


i/>  =  U-  ST 


where  U  and  S  respectively  denote  the  internal  energy  and  entropy.  As  detailed  in  [114],  the  combi¬ 
nation  of  mean  field  theory  and  classical  statistical  mechanics  arguments  yields  the  Helmholtz  energy 
relation 


4(P,T) 


T0iV 

2V 


[1  -  ( P/Ps )2]  + 


TkN 
2 VPS 


EhPs 

2 


[1  -  (P/Ps)2\  + 


Pln(^^)+Ps\n(l-(P/Ps)2) 

Pln(^^j+Psln(l-(P/Ps)2) 


(4) 


Here  Eh  =  y p2-  is  a  bias  field,  Tc  =  -y  is  the  Curie  temperature  (k  is  Boltzmann’s  constant),  and 
Ps  =  ^  denotes  the  saturation  polarization.  It  is  illustrated  in  [114]  that  i/j  yields  a  double  well 
potential  for  T  <TC  and  a  single  well  potential  for  T  >  Tc  in  accordance  with  the  transition  from  a 
ferroelectric  phase  to  a  paraelectric  phase. 
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Figure  3.  (a)  Helmholtz  energy  ifj  and  Gibbs  energy  G  for  increasing  fields  E.  (b)  Dependence 
of  the  local  average  polarization  P  on  the  field  E  at  the  lattice  level  in  the  presence  of  thermal 
activation. 

To  simplify  subsequent  implementation  for  fixed  temperature  regimes  T  <  Tc,  we  also  employ 
the  piecewise  quadratic  Helmholtz  relation 

'  ^v(p  +  Pr)2  ,P<~Pi 

i/;(P)  =  <  lv(P  -  Pr )2  ,  p>  Pi  (5) 

_  Hpi  ~  pr)  (w  -  Pr)  >  lpl  <  pi 

obtained  through  Taylor  expansion  of  (4)  in  neighborhoods  of  the  equilibria.  The  parameters  Pr  and 
Pj  respectively  denote  the  local  polarization  value  at  which  a  local  minimum  of  ^  occurs  (see  (1)) 
and  the  inflection  point.  From  a  physical  perspective,  Pr  and  r/  represent  the  local  renranence  value 
and  reciprocal  of  the  slope  at  saturation  for  the  hysteron  depicted  in  Figure  3b.  Finally,  we  note 
that  the  local  coercive  field  Ec  is  related  to  Pr,  Pj  and  rj  through  the  relation 

Ec  =  V(Pr-Pj).  (6) 

From  the  potential  energy  relation  Ue  =  —  p  •  E,  it  follows  from  (2)  that  an  expression  for  the 
Gibbs  energy  is 

G(P,  E,  T)  =  ip(P,  T)  —  EP  (7) 

for  uniaxial  lattice  structures.  The  behavior  of  G  for  an  increasing  field  and  the  resulting  local 
average  polarization  P  are  depicted  in  Figure  3. 

2.1.2  Local  Average  Polarization 

To  quantify  P  it  is  necessary  to  determine  the  fractions  of  dipoles  x+  and  X-  respectively  having 

positive  and  negative  orientations,  the  likelihoods  p. \ _ and  p _ of  switching  from  positive  to  negative, 

and  conversely,  and  the  expected  polarization  values  (P+)  and  (P~)  associated  with  positive  and 
negative  dipoles.  In  the  absence  of  thermal  activation,  these  quantities  follow  directly  from  the 
necessary  condition  (3)  which  forms  a  basis  for  subsequent  asymptotic  relations  for  P.  To  incorporate 
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measured  relaxation  phenomena,  it  is  necessary  to  incorporate  a  balance  of  the  thermal  and  Gibbs 
energies  through  the  Boltznran  relation 


H(G)  =  Ce~GV'kT. 


(8) 


The  constant  C  is  chosen  to  guarantee  integration  to  unity  over  all  admissible  inputs  and  k  again 
denotes  Boltzmann’s  constant.  Physically,  the  probability  p(G)  of  achieving  an  energy  level  G  is 
increased  for  large  values  of  the  relative  thermal  energy  kT /V  and  diminished  for  small  kT /V .  Hence 
more  dipoles  will  have  the  energy  required  to  jump  from  one  stable  equilibrium  to  the  other  when 
kT /V  is  large  than  when  kT /V  is  small  or  in  the  limiting  case  when  kT /V  — >  0. 

As  detailed  in  [114],  the  likelihoods  p y _ and  p _ are  given  by 


^-GiE^i^^V/kT 


0—G(E,—PI(T),T)V/kT 


P+-  = 


T(T) 


L 


e-G(E,P,T)V/kTdp 


P-+  = 


T(T)  r-Pi(T) 


(9) 


Pi(T) 


e-G(E,P,T)V/kTdp 


where  the  formulation  in  terms  of  the  inflection  points  ±P/  rather  than  the  equilibrium  value  Pq 
can  be  argued  either  from  energy  principles  or  asymptotic  analysis.  The  relaxation  time  T  for  the 
material  is  related  to  the  temperature  and  reference  volume  V  of  mass  v  by  the  relation 


T(T)  =  V  ~kE~  <10) 

which  can  be  interpreted  as  the  reciprocal  of  the  frequency  at  which  dipoles  attempt  jumps.  By 
integrating  the  product  Pp(G)  over  all  admissible  polarization  values  and  evaluating  the  integration 
constant  C,  it  follows  that  the  average  polarization  values  can  be  expressed  as 


/»oo 

/  p  e~G{E  ,P,T)V /  kT  dp 

( P+)  =  —r bo -  ,  (P-)  = 

/  e-G(E,P,T)V/kTdp 
JPt 


—  P, 


pe-G(E,P)V/kTdp 


>  — oo 

r-Pj 


-G(E,P)V/kTdp 


(11) 


The  evolution  of  dipole  fractions  is  quantified  by  the  differential  equations 


X+  =  -p+-X+  +P-+X- 
X-  =  —p _ yX-  +  P- 1 _ X+ 


(12) 


which  can  be  simplified  to 

X+  =  -P+-X+  +P-+(1  -  x+)  (13) 

through  the  identity  x+  +  X-  =  1. 

With  the  dipole  fractions,  transition  likelihoods,  and  average  polarization  values  thus  defined, 
the  average  local  polarization  at  the  lattice  level  is  given  by 


P  =  x+  (P+)  +  am  (P_) .  (14) 

The  formulation  (14)  incorporates  thermal  relaxation  through  the  relation  (10)  as  well  as  certain 
temperature  dependencies  if  the  Helmholtz  relation  (4)  is  employed.  While  the  use  of  (4)  has  been 
illustrated  for  incorporating  the  temperature  dependence  associated  with  certain  relaxor  ferroelectric 
compounds,  more  general  temperature  dependencies  can  be  incorporated  by  modifying  the  piecewise 
quadratic  relation  (5)  in  the  manner  illustrated  in  Section  2.3  for  shape  memory  alloys. 


2.1.3  Thermal  Evolution 


To  quantify  changes  in  temperature  due  to  convection  and  conduction  to  surrounding  media, 
Joule  heating,  and  dipole  switching,  we  employ  a  balance  of  energy  to  obtain  the  evolution  equation 

cMT{t)  =  -hcn[T  -  TE{t )]  —  -J~[T  ~  TE(t)}  +  J{t)  -  [h+x+  +  h-±-\.  (15) 

Here  c,Ai,  hc,  H,  A  and  i  respectively  denote  the  specific  heat  for  PZT,  the  mass  of  the  actuator,  a 
heat  transfer  coefficient,  the  surface  area  of  the  PZT,  the  thermal  conductivity  of  the  surrounding 
medium,  and  the  interval  over  which  conduction  occurs  [47] .  The  first  term  on  the  right  hand  side  of 
(15)  quantifies  heat  exchange  due  to  convection  whereas  the  second  term  incorporates  potential  heat 
loss  or  sources  due  to  conduction.  The  term  J{t )  characterizes  temperature  changes  due  to  Joule 
heating.  For  certain  operating  regimes,  this  can  be  quantified  by  the  relation 

J(t)  =  ifh!^  (16) 

where  pe,h  and  A  respectively  denote  the  average  electrical  resistivity,  thickness  of  the  PZT,  and 
cross-sectional  area  of  the  PZT  [100].  To  incorporate  additional  geometric,  electromagnetic,  or 
frequency-dependent  effects,  one  can  quantify  J(t)  either  through  experimental  measurements  or 
more  comprehensive  models.  We  note  that  the  incorporation  of  Joule  heating  mechanisms  becomes 
increasingly  important  in  applications  requiring  high  drive  frequencies. 

The  final  component  of  (15)  quantifies  heat  transduction  due  to  dipole  switching  so  that  h+  and 
h-  are  analogous  to  the  specific  enthalpies  in  the  corresponding  SMA  relation  (51).  This  relation  also 
becomes  increasingly  significant  for  high  frequency  transduction.  In  general,  validation  experiments 
will  be  required  to  establish  regimes  when  this  latter  contribution  should  be  retained  as  well  as 
operating  conditions  where  it  can  be  considered  negligible.  Finally,  we  point  out  that  c  and  pe  are 
considered  as  averages  and  the  specification  of  phase-dependent  components  to  these  coefficients  may 
be  required  when  quantifying  the  temperature  changes  which  occur  during  phase  transitions. 


2.1.4  Asymptotic  Polarization  Relations  in  the  Absence  of  Thermal  Activation 


For  applications  in  which  relaxation  times  are  considered  negligible,  or  do  not  constitute  a  sig¬ 
nificant  role  in  transducer  behavior,  the  local  average  polarization  expression  can  be  simplied  signif¬ 
icantly  by  considering  the  limiting  case  in  which  thermal  fluctuations  about  equilibrium  values  are 
eliminated.  As  detailed  in  [114],  the  limiting  models  can  be  derived  either  through  formal  application 
of  the  necessary  condition  (3)  or  rigorous  consideration  of  the  limits  kT /V  — »  0  in  the  definition  of 

the  likelihoods  p. \ _ ,  p _ (.  and  average  polarization  values  (P+) ,  (P_).  To  simplify  the  discussion,  we 

summarize  the  formal  arguments  here. 

For  the  Helmholtz  relation  (4),  enforcement  of  the  necessary  condition  (3)  with  e  =  P  and  ip  =  E 
yields  the  equilibrium  relation 

P{E)  =  Pstanh  >  (17) 


where 


a(T) 


EhT 

Ps 


for  the  average  local  polarization.  We  note  that  (17)  is  exactly  the  Ising  relation  employed  in  [108] 
when  quantifying  the  anhysteretic  polarization  and  is  equivalent  through  first-order  terms  with  the 
Langevin  relations  employed  in  various  ferroelectric  and  ferromagnetic  models.  Translates  of  the 
form  P  =  Pstanh(P  ±  Ec )  were  also  employed  by  Zhang  and  Rogers  [123],  and  ridge  functions 
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of  the  form  r(x)  =  tanh(x')  are  appropriate  for  generalized  Preisach,  or  Krasnosel’skii-Pokrovskii 
characterizations  [8,  9,  60]. 

The  relation  resulting  from  (5)  is  slightly  more  complicated  due  to  the  piecewise  nature  of  the 
definition.  Enforcement  of  the  necessary  condition  (3)  yields  the  general  relation 

P=-E  +  A  PR  (18) 

V 

where  A  =  1  for  positively  oriented  dipoles  and  A  =  —  1  for  negative  orientations.  To  specify  P  more 
rigorously  in  terms  of  initial  dipole  orientations 

'  f  ~Pr  ,  E( 0)  <  -Ec 

[P(E;Ec,0m=  £  ,-Ec<E(0)<Ec  (19) 

,f  +  Pr  ,  E{ 0)  >  Ec 

and  transition  times 

r{t)  =  {t  G  (0,  Tf]  |  E(t)  =  -Ec  or  E(t)  =  Ec},  (20) 

we  employ  the  Preisach  notation 

'  \P(E-,Ec,Om  *r(t)  =  0 

[P(E;Ec,£)](t)  =  <  f  -  Pr  ,  r(t)  /  0  and  £(maxr(f))  =  -Ec  (21) 

,  f  +  Pr  ,  r(t)  /  0  and  E(maxr(t))  =  Ec 

(e.g.,  see  [8,  60,  112]).  The  dependence  of  the  local  average  polarization  P  on  the  local  coercive  field 

Ec  is  indicated  as  a  prelude  to  the  stochastic  extensions  to  the  theory  developed  in  Section  3. 

2.2  Magnetostrictive  Materials 

Field-induced  magnetization  changes  in  ferromagnetic  compounds,  and  accompanying  hysteresis, 
are  due  primarily  to  a  combination  of  two  mechanisms:  rotation  of  moments  and  losses  associated 
with  domain  wall  motion  [19,  27,  56].  For  highly  anisotropic  materials,  characterization  of  consti¬ 
tutive  nonlinearities  and  hysteresis  necessarily  requires  incorporation  of  the  anisotropy  energy  when 
computing  the  energetic  response  of  the  material.  For  several  classes  of  operating  conditions  and 
materials,  however,  models  based  on  the  restriction  of  spins  or  magnetic  moments  to  two  possible  ori¬ 
entations,  denoted  Sj  =  ±1,  can  be  developed.  From  a  classical  perspective,  this  assumption  is  valid 
for  materials  which  exhibit  uniaxial  crystalline  anisotropies  or  transducers  in  which  uniaxial  stresses 
dominate  the  crystalline  structure.  This  latter  category  includes  Terfenol-D  transducers  in  which 
prestress  mechanisms  are  employed  to  maintain  the  rod  in  compression  and  to  optimize  magnetonre- 
chanical  coupling  [29,  30].  From  a  quantum  perspective,  the  assumption  of  diametrically  opposed 
spins  is  commensurate  with  the  observation  that  two  allowable  spin  configurations  are  orientations 
parallel  and  opposite  to  an  applied  field. 

2.2.1  Helmholtz  and  Gibbs  Energies 

To  construct  a  mesoscopic  Helmholtz  energy  relation  for  ferromagnetic  materials,  we  consider  a 
uniform  and  homogeneous  lattice  volume  V  of  mass  v  comprised  of  N  cells,  each  of  which  is  assumed 
to  contain  one  spin,  or  magnetic  moment,  m  having  the  orientation  Sj  =  ±1.  As  in  the  ferroelectric 
model,  we  let  $0  denote  the  energy  required  to  reorient  a  single  moment  in  a  completely  ordered 
lattice  -  the  relation  between  $0  and  the  magnetic  exchange  integral  J  is  detailed  in  [105]  -  and  make 
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the  assumption  that  only  adjacent  moments  interact  in  accordance  with  mean  field  approximation 
to  the  Ising  model.  As  detailed  in  [105],  formulation  of  the  Helmholtz  energy  ip  =  U  —  ST  in  this 
context  yields  the  relation 


HhMs  n  (1/ri1/r\ 21  ,  HhT 


M  In 


f  M  +  Ms\ 
V  Ms  —  M  ) 


+  MS  ln(l  -  (M/Ms)2) 


(22) 


where  Ms  =  Tc  =  and  Hy  =  respectively  denote  the  technical  saturation  magnetization, 
Curie  temperature,  and  bias  field.  The  initial  assumption  that  4>o  is  constant  implies  that  Hf l  is 
constant  for  a  uniform,  homogeneous  lattice.  This  assumption  is  relaxed  in  Section  3  where  we 
consider  statistically  distributed  values  of  H y  to  accommodate  material  nonhomogeneities  and  the 
effects  of  polycrystallinity.  Finally,  we  note  that  ip  yields  a  single  well  potential  for  T  >  Tc  and  a 
double  well  potential  in  accordance  with  the  transition  from  paramagnetic  to  ferromagnetic  phases. 

In  [105],  it  is  illustrated  that  Taylor  expansion  about  the  three  equilibria  can  be  employed  to 
construct  the  piecewise  quadratic  Helmholtz  relation 


ip(M) 


'  \q{M  +  MR )2 
<  \v{M-Mr)2 


,  M  <  -Mr 
,  M  >  Mr 

,  \M\  <  Mr 


(23) 


which  facilitates  implementation  in  fixed  temperature  regimes.  Here  Mj  and  Mr  denote  the  inflection 
point  and  minimum  of  ip  at  which  local  remanence  occurs  -  see  Figure  3  for  a  depiction  of  analogous 
behavior  for  ferroelectric  materials. 

For  variable  temperature  applications,  one  can  employ  either  (22)  or  a  modification  of  (23)  in 
which  temperature  dependence  is  incorporated  in  a  manner  analogous  to  that  employed  in  Section  2.3 
for  shape  memory  alloys.  The  former  option  automatically  includes  the  ferromagnetic  to  paramag¬ 
netic  phase  transition  whereas  the  latter  offers  advantages  from  the  perspective  of  implementation. 

Since  the  magnetostatic  energy  is  given  by  £  =  —  Hom  •  H,  where  ho  denotes  the  magnetic 
permeability,  the  uniaxial  Gibbs  energy  can  be  specified  either  as 


G(M,  H ,  T)  =  iP(M,  T)  -  hq MH 


or 

G(M,  H,  T)  =  ip(M,  T)  -  MH  (24) 

by  incorporating  ho  in  V;-  We  employ  (24)  due  to  its  commonality  with  the  general  Gibbs  relation  (2) 
which  facilitates  the  development  of  a  unified  modeling  framework.  A  comparison  of  the  Helmholtz 
relations  (4)  and  (22),  the  piecewise  quadratic  relations  (5)  and  (23),  and  the  Gibbs  energy  relations 
(7)  and  (24)  illustrates  the  commonality  in  the  formulations  for  ferroelectric  and  ferromagnetic 
compounds. 

2.2.2  Average  Local  Magnetization 

The  development  of  local  average  magnetization  relations  is  analogous  to  that  for  ferroelectric 
materials  and  details  are  provided  in  [105,  106].  For  operating  regimes  in  which  relaxation  effects  are 
considered  negligible,  limiting  values  of  the  local  average  magnetization  M  can  be  formally  obtained 
from  the  necessary  condition  (3).  For  the  Helmholtz  relation  (22)  derived  from  statistical  mechanics 
principles,  the  condition  H  =  yields 

M  =  ATTanh  (25) 
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TT  JJ  ( 'J A 

where  a  =  jp  and  a(T )  =  ^  '  .  The  relation  (25)  is  an  Ising  model  whose  input  is  the  effective 

field 

He  =  H  +  aM.  (26) 

This  relation  was  employed  as  the  anhysteretic  component  of  the  unified  models  developed  in  [109]. 
Furthermore,  it  is  illustrated  in  [108,  109]  that  if  one  relaxes  the  constraint  that  moments  have  only 
the  orientations  s*  =  ±1  and  considers  uniformly  distributed  moments,  one  obtains  the  Langevin  rela¬ 
tion  M  =  C(He)  =  Ms[coth(He/a )  —  a/He ]  which  agrees  with  the  Ising  relation  M  =  Mstanh(i7e/a) 
through  first  order  terms  -  see  [19,  56]  for  a  derivation  of  the  Langevin  equation  in  the  context  of 
magnetic  materials.  The  Langevin  model  M  =  MsC(He),  with  He  specified  by  (26),  is  employed 
when  quantifying  the  anhysteretic  magnetization  in  the  domain  wall  theory  of  Jiles  and  Atherton 
[57,  58]  as  well  as  the  transducer  models  based  on  that  theory  [17,  29,  30]. 

The  magnetization  relation  resulting  from  (24)  in  the  absence  of  thermal  activation  is  elementary 
in  the  sense  that  it  is  piecewise  linear  but  is  complicated  by  the  fact  that  a  history  of  dipole  switches 
must  be  maintained  to  ascertain  which  branch  of  the  hysteron  is  active.  This  can  be  accomplished 
via  the  Preisach  notation 


[M(H-Hc,Om=  { 


(  [M(H;Hc,Om  ,  T(t)  =  0 

—  Mr  ,  r(t)  /  0  and  H (max  r(t))  =  —  Hc 

+  Mr  ,  r(t)  /  0  and  H (max  r(t))  =  Hc 


(27) 


where 


(  H 


[M(H;HC,  £)](0)  =  { 


f  -  Mr  ,  H( 0)  <  -Hc 


£ 

H 


,  -Hc  <  H( 0)  <  Hc 


(28) 


l  f  +  Mr  ,  H( 0)  >  He 


denotes  initial  moment  orientations  and  the  switching  times  r(t)  are  specified  by  (20)  with  switching 
occurring  at  the  local  coercive  field  values  Hc  and  —Hc.  Heuristically,  the  relation  (27)  has  the  form 


M  =  -H  +  A  Mr 
V 


(29) 


where  A  =  ±1  is  used  to  delineate  the  upper  and  lower  branches  of  the  hysteron. 

To  incorporate  thermal  activation  and  relaxation  mechanisms  -  e.g.,  when  modeling  magnetic 
after-effects  of  the  type  discussed  in  Chapter  20  of  [19]  -  it  is  necessary  to  quantify  the  evolution  of 
the  moments  fractions  and  X-  having  positive  and  negative  spins.  The  notation  and  development 
is  analogous  to  that  presented  in  Section  2.1.2  for  ferroelectric  materials  and  details  regarding  the 
development  of  ferromagnetic  materials  can  be  found  in  [105,  106]. 

Since  the  probability  of  achieving  an  energy  state  G  is  given  by  (8),  the  transition  likelihoods 
p. ) _ and  p _ are  given  by 


P+-  = 


T(T) 


roo 
J  M, 


0—G(H,Mi,T)V/kT 
e~G(H,M,T)V/kT  dM 


P-+  = 


T(T) 


e 

r—Mj 


-G(H-M!,T)V/kT 


-G(H,M,T)V/kTdM 


(30) 


where  the  relaxation  time  T  is  specified  by  (10).  The  expected  magnetization  values  (M+)  and 
(M_),  due  to  positively  and  negatively  oriented  moments,  are  specified  by 


(M+)  = 


poo 
JM t 


M  e~G(H,M,T)V/kT  dM 


»- Mj 


M  e~G(H,M)V/kT  dM 


poo 

JMj 


-G(H,M,T)V/kT  dM 


(M~)  = 


-G(H,M)V/kT  dM 


(31) 
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which  are  analogous  to  the  ferroelectric  relations  (11).  Finally  the  average  local  magnetization  A f, 
which  incorporates  thermal  relaxation  mechanisms,  is  given  by 

M  =  x+  (M+)  +  X-  (M_)  (32) 

where  the  moment  fraction  is  specified  by  the  differential  equation  (13)  and  x_  =  1  —  x+.  In 
[105],  it  is  illustrated  that  the  relations  (27)  for  materials  with  negligible  relaxation  characteristics 
can  also  be  rigorously  derived  from  (30),  (31)  and  (32)  by  considering  the  limiting  case  kT /V  — >  0. 

2.2.3  Thermal  Evolution 

Temperature  changes  in  magnetostrictive  transducer  materials  are  due  to  several  mechanisms 
including  ohmic  heating  in  the  solenoid,  which  is  transmitted  to  the  Terfenol-D  rod  through  con¬ 
duction,  Joule  heating  due  to  eddy  currents,  and  potential  internal  heating  due  to  the  reorientation 
of  moments.  To  accommodate  future  designs,  we  also  incorporate  possible  convection  mechanisms 
analogous  to  those  considered  for  piezoceramic  and  shape  memory  alloy  transducers. 

To  model  the  temperature  transition  in  the  magnetic  material,  we  employ  the  notation  of  Sec¬ 
tion  2.1.3  and  let  c,  hc,  AA,  fi,  A,  l  and  TA(t)  respectively  denote  the  specific  heat,  a  convection  co¬ 
efficient,  the  mass  of  the  Terfenol-D  rod,  the  surface  area  of  the  material,  the  thermal  conductivity 
of  the  solenoid,  the  path  length  of  conduction,  and  the  time  varying  temperature  of  the  adjacent 
environment.  An  energy  balance  then  yields  the  differential  equation 

cMT(t)  =  -n  [hc  +  A /£}  [T  -  Te{T)}  +  J(t)  -  [h+x+  +  h-i-]  (33) 

governing  the  temperature  evolution.  A  comparison  of  (33)  with  its  ferroelectric  counterpart  (15) 
reveals  that  the  general  temperature  relations  are  identical. 

To  quantify  the  conductive  and  convective  contributions,  it  is  necessary  to  specify  Tg(t)  either 
experimentally  or  through  additional  models.  For  example,  a  fully  coupled  thermal  model  for  a 
Terfenol-D  transducer  would  require  the  measurement  or  characterization  of  ohmic  heating  in  the 
solenoid  and  heat  transmission  through  the  device  to  specify  T^(t).  Similarly,  the  Joule  component 
J(t)  can  be  determined  either  empirically  or  though  eddy  current  relations.  Finally,  the  last  term  in 
(33)  quantifies  potential  changes  in  temperature  due  to  the  reorientation  of  moments.  The  relevance 
of  this  term  for  specific  operating  conditions  should  be  established  through  validation  experiments. 

2.3  Shape  Memory  Alloys 

Shape  memory  alloys  exhibit  a  number  of  unique  features  which  are  being  targeted  for  present  and 
projected  applications  [69,  99,  101].  Following  a  plastic  deformation  at  low  temperatures,  the  materi¬ 
als  will  recover  their  original  shape  upon  the  application  of  heat  -  this  constitutes  the  shape  memory 
effect  which  gives  the  materials  their  name.  At  higher  temperatures,  reversible  deformations  up  to 
10%  can  be  obtained  under  nearly  constant  loads  which  gives  the  materials  superelastic  capabilities. 
These  deformation  properties  are  commonly  termed  quasiplastic  at  low  temperatures  and  pseudoe¬ 
lastic  at  higher  temperatures  with  the  latter  illustrated  by  the  SMA  data  plotted  in  Figure  lc.  All 
of  these  characteristics  are  due  to  temperature  or  load  induced  phase  transformations  at  the  lattice 
level  and  it  is  at  this  scale  that  we  initiate  model  development.  Most  initial  applications  utilizing 
SMA  compounds  were  focused  on  polycrystalline  NiTi  wires  whereas  NiTi  films  are  under  present 
investigation  to  support  MEMs  applications  and  to  improve  response  times  through  increased  surface 
to  volume  ratios.  For  uniaxial  loads,  uniaxial  models  suffice  for  both  regimes  and  are  the  focus  of 
this  discussion.  Hence  we  focus  on  characterizing  the  evolution  of  austenite  A  and  the  martensite 
twins  Al+  and  M_  in  SMA  lattice  layers. 


13 


Following  the  theory  developed  in  [69,  79,  93,  98,  99],  we  treat  a  lattice  volume  V  of  mass  v  as 
the  fundamental  element  in  the  model  and  let  x+(f)  and  x-(t)  respectively  denote  the  volume 

fraction  of  A,  M+  and  M~  layers  in  the  SMA.  The  phase  fractions  constitute  internal  variables  which 
necessarily  satisfy  the  conservation  relation 

xa  +  x+  +  X-  =  1  (34) 

over  all  time.  To  consolidate  notation,  we  let  a  generically  refer  to  the  A,  M+  and  M~  variants  so 
that  the  relation  (34)  can  be  reformulated  as 

^2,xa{t)  =  i. 


2.3.1  Helmholtz  and  Gibbs  Energies 

The  formulation  of  Helmholtz  and  Gibbs  energy  relations  is  common  to  a  number  of  SMA  theories 
(e.g.,  see  [53,  69,  99,  101]  and  the  references  cited  therein),  and  differences  between  the  theories  arise 
in  the  manner  through  which  energy  relations  are  defined  and  used  to  construct  macroscopic  models. 
We  employ  the  theory  of  [69,  99]  to  construct  C\  energy  relations  based  on  piecewise  quadratic 
functions  which  characterize  stable  equilibria  corresponding  to  the  A,  M+  and  M -  phases. 

As  illustrated  in  Figure  4,  the  fundamental  order  parameter  is  the  shear  strain  e  which  has  the 
value  s  =  0  for  austenite  and  the  equilibrium  value  £  =  £t  for  martensite  in  a  stress-free  state.  The 
linear  moduli  for  martensite  and  austenite  are  respectively  denoted  by  Ym  and  Ya ■  To  model  the 
transition  from  stability  of  martensite  variants  at  low  temperatures  to  stability  of  the  austenite  phase 
at  high  temperatures,  we  employ  the  C\  Helmholtz  energy  relation 

£  <  —£m(T) 

—£m{T)  <  £  <  —£a{T) 

kl  <  sa(T)  (35) 

£a{T)  <  £  <  £m(T) 

£  A  £m(T)  ■ 

and  their  negatives,  delineate  the  transi¬ 
tion  from  convex  regions,  which  represent  stable  austenite  and  martensite  phases,  to  concave  regions 

8  =  £t 


M“  A  M+ 

Figure  4.  Lattice  element  exhibiting  the  martensite  M~,M+  and  austenite  A  equilibrium  configu¬ 
rations. 


iP(£,t)  =  ! 


im  ,  ,  o 

~Y^  +  £t) 


E0(T) 


Ya_2 


(£  +  £0(T))2  +  MT) 


~fe2  +  A  (3(T) 


Eo(T) 


(■ s-£0(T))2  +  MT ) 


Ym  ,  s2 


The  temperature-dependent  inflection  points  £m{T),£a{T), 
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Figure  5.  (a)  Piecewise  quadratic  Helmholtz  energy  (35)  for  fixed  temperature  T ;  dashed  segments 
represent  concave,  unstable  regions,  (b)  Helmholtz  energy  as  a  function  of  T. 

representing  unstable  states.  The  maxima  of  the  concave  parabolae  occur  at  the  temperature- 
dependent  points  (±eo(T) ,  ipo(T)) ,  and  the  parameter  Eq(T)  is  chosen  to  ensure  C\  continuity.  As 
depicted  in  Figure  5,  the  austenite  minimum  has  the  height 

A (5{T)  =  pA(T)  -  pM(T)  (36) 

where 

(3a(T)  =  ca(T  -  Tr)  +  uR  -  Tsa  (37) 

represent  chemical  (nonelastic)  free  energies  [69,  99,  101].  Here  ca,ua  and  Tr  respectively  denote 
specific  heat  capacities  multiplied  by  the  material  density,  phase-dependent  internal  energy  constants, 
and  the  temperature  of  the  reference  state  from  which  energies  are  computed.  Furthermore, 

Sa  =  ca  In (T/Tr)  +  r]R  (38) 

are  specific  entropies  and  r]Q  are  phase-dependent  entropy  constants. 

The  Helmholtz  relation  (35)  quantifies  the  energy  for  a  layer  in  the  absence  of  applied  loads. 
To  incorporate  distortions  in  the  energy  landscape  due  to  applied  stresses  a,  we  employ  the  Gibbs 
relation 

G(e,  a,  T)  =  -0(e,  T)  —  ecr  (39) 

resulting  from  (3).  In  Figure  6a,  we  illustrate  the  Gibbs  energy  for  a  lattice  volume  V  at  a  fixed 
temperature  T  =  Ta  chosen  so  that  the  material  exhibits  the  austenite  phase  for  a  =  0.  As  a 
is  increased,  the  landscape  distorts  until  the  critical  stress  a  a  at  which  point  the  stable  austenite 
equilibrium  ceases  to  be  a  local  minimum  and  M+  becomes  the  stable  phase.  It  will  remain  such  until 
the  stress  is  decreased  to  a  second  critical  value  &M-  Here  the  local  martensite  minimum  disappears 
and  the  material  returns  to  the  austenite  phase. 

2.3.2  Local  Average  Stress-Strain  Relations 

To  quantify  the  pseudoelastic  load  deformation  behavior  which  accompanies  the  increase  in  stress 
through  the  critical  value  a  a  and  subsequent  decrease  through  the  second  critical  value  ctm>  one  can 
quantify  the  evolution  of  associated  phase  fractions  in  the  absence  or  presence  of  thermal  activation 
in  a  manner  analogous  to  that  summarized  previously  for  ferroelectric  and  ferromagnetic  compounds. 

For  reference  volumes  V  in  which  thermal  activation  is  negligible,  the  local  average  stress- strain 
behavior  is  quantified  solely  by  necessary  conditions  associated  with  the  Gibbs  energy  G  given  by  (39). 
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Figure  6.  (a)  Gibbs  energy  (39)  for  a  fixed  temperature  T  in  the  austenite  range  with  o  =  0, 
<7  =  a  a  and  a  =  a  m-  (b)  Local  stress-strain  relation  produced  by  the  equilibrium  condition  (3)  for 
fixed  temperature  T  and  negligible  relative  thermal  energy  kT/V.  (c)  Local  stress-strain  relation  in 
the  presence  of  relative  thermal  energy  kT/V  for  fixed  T. 


From  (2),  it  follows  that  the  conditions 


dG 


=  0 


d2G 


>  0 


de  ”  ’  de2 

must  hold  at  stable  equilibria  which  yields  the  critical  stress  values 

aA(T)  =  YAeA(T) 


(40) 


(41) 


&m{T)  =  Ym[sm{T)  -  £t\. 

For  material  characterization,  it  is  advantageous  to  define  the  temperature-dependent  parameters 

S(T)  =  aA(T )  -  aM(T)  (42) 


which,  it  will  be  later  illustrated,  can  be  estimated  from  measured  data.  Due  to  the  quadratic 
definition  of  the  Gibbs  energy,  the  local  stress-strain  behavior  is  linear  in  the  absence  of  thermal 
activation,  as  depicted  in  Figure  6b.  To  incorporate  the  switching  history,  we  employ  the  Preisach 
notation 


[e(M,  £)](*) 


[e(cr;<5,£)](0)  ,  T(t)  =  0 

<  +  £t  ,  r(t)  /  0  and  ff(maxr(t))  =  aA 

y-  ,  r(t)  /  0  and  ff(maxr(t))  =  ctm 


(43) 


to  quantify  the  local  average  strains  due  to  the  positive  applied  loads  a  (similar  expressions  hold 
for  compressive  loads).  Analogous  to  the  ferroelectric  definitions  (19)  and  (20)  and  ferromagnetic 
definition  (28), 


Ya 

£ 


l  Ym 


+  £t 


[e(<r;<J,C)](0)  =  < 
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o-(O)  <  aM 
aM  <  cr(0)  <  a  a 
o-(O)  >  aA 


(44) 


and 


(45) 


r(t)  =  {t  <E  (0,  Tf]  I  a(t)  =  aM  or  a(t)  =  aA}, 

respectively  denote  the  initial  phase  configuration  and  switching  times. 

To  quantify  the  local  average  strains  e  for  reference  volumes  in  which  thermal  activation  is  signif¬ 
icant,  it  is  necessary  to  balance  the  Gibbs  energy  G  with  the  relative  thermal  energy  kT/V  through 
the  Boltzmann  relation  (8).  Physically,  large  values  of  kT/V  mean  that  with  higher  probability, 
layers  will  achieve  the  energy  required  to  exit  a  local  minimum  before  the  stress  values  a  a  or  ctm 
are  reached.  This  produces  the  gradual  transitions  and  decreased  transition  stress  values  depicted 
in  Figure  6c. 

As  detailed  in  [69,  93,  99],  the  likelihood  pA±  that  austenite  will  transform  to  M±  and  the 
likelihoods  p±  that  M±  will  transform  either  to  austenite  or  the  other  martensite  variant  are  given  by 

e-G(±6A,a,T)V/kT 
r  e-G(e,a,T)V/kTde 

JXA 

(46) 

e-G(±sM,<r,T)V/kT 

f  e-G{e,a,T)V/kTd£ 

JXM± 

where  xa{T)  =  (-£A(T),ea(T)),xm+(T)  =  (eM(T),  oo)  and  Xm~{T)  =  (-00 ,-eM(T))  respectively 
denote  regions  over  which  austenite,  M+  and  M~  are  stable.  Moreover,  the  expected  strains  due  to 
austenite  and  M±  variants  are  given  by 


P±(v,T)  = 


I  kT 
2ttuV 2/3 


pA±{a,T )  = 


I  kT 
2-KuV2/‘i 


<£-)  = 


e//(G(e,  <7,  T))  de 


{e+)  = 


sn(G(e,  a,  T))  de 


(eA)  =  /  ep(G(e,a,T))de 

JXA 

where  p  is  given  in  (8)  and  the  Gibbs  energy  is  specified  in  (39).  The  evolution  of  phase  fractions  is 
governed  by  the  rate  laws 

X-(t)  =  pA-XA(t)  P—X—  (t) 

x+(t)  =  pA+xA{t)  -  p+x+(t)  (48) 

xA{t)  =  P-X-(t )  -  pA-XA(t )  +  p+x+{t)  -  pA+xA{t ) 
which  can  be  reduced  to 

x-(t)  =  -{p-  +  pA—)x—(t)  -pA-X+(t)  +pA- 

(49) 

x+(t)  =  ~(p+  +pA+)x+(t )  -pA+X-(t)  +pA+ 

through  the  conservation  relation  (34).  To  compute  the  local  average  strain  for  the  reference  volume, 
we  make  the  assumption  that  thermal  strains  are  small  compared  to  mechanical  strains  and  retain 
only  the  latter  component.  As  discussed  in  [69],  this  assumption  is  valid  for  bulk  materials  but  may 
need  to  be  modified  for  SMA  thin  films.  Under  this  assumption,  the  local  average  strains,  in  regimes 
for  which  thermal  activation  or  relaxation  mechanisms  are  significant,  are  given  by 

e  =  (e-)x-  +  (e+)x+  +  {eA)xA  (50) 

where  x-,x+  are  specified  by  (49),  xA  =  1  —  x+  —  X-,  and  (e_) ,  (e+)  and  (eA)  are  defined  in  (47). 


17 


2.3.3  Thermal  Evolution 

The  quantification  of  thermal  processes  is  analogous  to  that  considered  in  (15)  for  ferroelectric 
compounds  and  (33)  for  ferromagnetic  materials,  and  includes  convective  and  conductive  mecha¬ 
nisms,  Joule  heating  and  heat  transduction  due  to  phase  transitions.  As  in  Sections  2.1.3  and  2.2.3, 
we  let  hc,  Te,  pea ,  ca,  0,  A,  At,  A  and  i  respectively  denote  the  convection  coefficient,  temperature 
of  the  surrounding  environment,  phase-dependent  electric  resistivity,  phase-dependent  specific  heat, 
SMA  surface  area,  SMA  cross-sectional  area,  the  mass  of  the  SMA,  the  thermal  conductivity  of  the 
surrounding  medium  and  the  length  of  any  conduction  paths.  An  energy  balance  then  yields  the 
differential  equation 


Mc(t)T(t)  =  -n[hc  +  \/l]  [T  -  Te(t)\  +  J(t)  -  ^  hQxa  (51) 

a 

which  quantified  the  temperature  evolution  (see  [69,  93,  98,  99]).  The  first  term  on  the  right  hand  side 
characterizes  heat  exchanged  with  the  surrounding  environment  through  convection  and  conduction 
in  a  manner  analogous  to  (15)  and  (33)  for  ferroelectric  and  ferromagnetic  materials.  For  an  input 
current  I(t),  the  relation 

=  (52) 

where  pe{t)  =  Yha  Paxa(t )  is  the  average  resistivity  per  unit  length,  quantifies  heat  generated  through 
Joule  heating.  The  final  term  in  (51)  accounts  for  heat  generated  or  lost  during  phase  transformations. 
As  discussed  in  [69],  the  specific  enthalpies  ha  have  the  form 

ha  =  9  a  +  Tsa 

where  ga  are  local  minima  of  the  Gibbs  relation  (39)  and  sa  denotes  the  specific  entropies  defined 
in  (38).  Finally,  the  average  specific  heat  is  given  by  c(t)  =  ^acaxa(t).  We  note  that  c  and  pe  are 
assumed  to  be  constant  in  the  ferroelectric  and  ferromagnetic  relations  whereas  their  time-dependence 
in  (51)  and  (52)  reflects  their  dependence  on  evolving  phase  distributions  in  shape  memory  alloys. 

2.4  Unified  Mesoscopic  Models 

In  Section  2.1,  2.2  and  2.3,  we  summarized  the  development  of  nonlinear  hysteresis  models  for  ferro¬ 
electric,  ferromagnetic  and  ferroelastic  compounds,  respectively.  We  discuss  here  the  unified  nature 
of  the  modeling  framework  and  indicate  the  degree  to  which  the  three  models  can  be  consolidated 
through  common  notation. 

As  noted  in  the  introduction  to  this  section,  the  polarization  P,  magnetization  M ,  and  strain 
e  can  all  be  represented  by  the  order  parameter  e  whereas  the  electric  field  E,  magnetic  field  H 
and  stress  a  can  be  respectively  represented  by  the  external  field  ip  which  is  thermodynamically 
conjugate  to  e.  The  definitions  of  subscripted  parameters  follow  the  conventions  established  in 
previous  subsections  (e.g.,  ej  =  Pi  or  Mi).  In  all  definitions,  the  Boltzmann  probability  n(G)  is 
specified  by  (8)  and  domains  of  integration  are  defined  by  %+  =  (ej,oo),x_  =  (— oo,  —  ei),XA  = 
(—£a,£a),Xm+  ={£m,  oo)  and  Xm-  =  (~ oo,  —  £m)- 

By  employing  the  order  parameter,  conjugate  fields,  and  appropriate  parameters,  the  magneti¬ 
zation  and  polarization  models  can  be  completely  unified  in  the  absence  of  phase  transitions.  The 
strain  model  for  ferroelastic  compounds  falls  within  the  same  modeling  framework  but  some  individ¬ 
ual  definitions  differ  slightly  due  to  the  inherent  phase  transitions.  While  notation  can  be  established 
to  formulate  all  three  models  in  terms  of  unified  definitions  and  expressions,  the  required  generality 
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obscures  rather  than  clarifies  the  framework  so  we  omit  it  here.  Finally,  we  note  that  the  uni¬ 
fied  mesoscopic  models  summarized  here  represent  only  the  fundamental  theory  and  a  number  of 
generalizations  to  this  theory  are  discussed  in  previous  subsections. 


Ferroelectric  and  Ferromagnetic  Materials 

Helmholtz  and  Gibbs  Energy  Relations: 


=  l-(e/e,)2]+|^ 


e,  —  e 


eln  [  6  +  6s  )  +  es  ln(l  -  (e/es)2^ 


■0(e) 


\r)(e  +  eR)2 
<  3??(e  -  eR)2 

\v(ei  ~  eR)  -  eR^ 


,  e  <  -e/ 
,  e  >  e7 

,  |e|  <  ei 


G(e,  tp,  T)  =  ip(e,T)  -  &p 


Mesoscopic  Model 

[e(<p;¥>c,£)](  0) 
r(t)  =  {te  (0,  Tf]  |  <p(t) 


,  r(t)  /  0  and  ^(maxr(t))  =  —fc 
,  r(i)  /  0  and  ^(maxr(t))  =  <£>c 

,  <^(0)  <  -<y9c 
,  -<^c  <  <^(0)  <  ipc 
,  ^(0)  >  Tc 

or  v?(t)  =  tfc } 


./Vo  Relaxation  Mechanisms: 

f  -  eR 

^-eR 
£ 


Mesoscopic  Model  -  with  Relaxation  Mechanisms: 

e  =  x+  (e+)  +  x_  (e_) 

(e±)  =  [  efi{G{e,  0,  T))de 
■'x± 

x+  =  -P+-X+ +  P-+{1  -  X+) 
x_  =  1  —  x+ 

P±  =  ^^y^(G(±e/,0,T)) 

Thermal  Evolution: 

cMT(t)  =  -n[hc  +  \/£\[T  -  TE{t )]  +  J(t)  -  /iaxa 
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Ferroelastic  Materials  —  Austenite  and  Martensite  Phases 

Helmholtz  and  Gibbs  Energy  Relations: 


V>(e,T) 


^(e  +  eT)2 

_^(e  +  eo(T))2  +  V>o(T) 

<  ^e2  +  A  (3(T) 

~^f^(e-e0(T))2  +  MT) 


,  e  <  -eM(T) 

,  -eM(T)  <  e  <  -eA(T) 
,  |e|  <  eA{T) 

,  eA(T)  <  e  <  eM{T) 

,  e  >  eM(T) 


G(e,ip,T)  =  ip(e,T)  -  ep 


Mesoscopic  Model  -  No  Relaxation  Mechanisms: 


[e(^;^c,?)](0) 

<p 

YA 


,  r(t)  =  0 

,  r(t)  /  0  and  <^(maxr(t))  =  pm 
,  r(t)  /  0  and  <p(maxr(t))  =  pA 


[e(¥>;^c,O](0) 


Ya 

Z 


,  p(0)  <  ipM 

,  pM  <  <fi(  0)  <  ip  A 

,  ¥>(0)  >  pA 


r(t)  =  {t  €  (0,  Tf\  |  p(t)  =  pM  or  p(t)  =  pA } 


Mesoscopic  Model  -  with  Relaxation  Mechanisms: 

e  =  x+  (e+)  +  x-  (e_)  +  xA  (eA) 

(e±>  =  [  en(G(e,  0,  T))de  , 

x-(t)  =  -{p-  +pA-)x-(t )  -pA-X+{t)  +PA- 
x+(t)  =  ~{p+  +  pA+)x+{t)  -pA+X-{t)  +pA+ 

XA  =  1  —  x_|_  —  ,T_ 

=  ^^y//(G(±eM,0,T))  ,  PA±  = 

Thermal  Evolution: 

McT{t)  =  -n[hc  +  \/£}[T  -  TE{t )]  +  J{t)  -  Y  Kxa 


(eA)  =  f  ep,(G(e,<j>,T))de 
Jxa 
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2.5  Macroscopic  Models  for  Homogeneous,  Single  Crystal  Compounds 

The  hysteresis  models  in  Sections  2.1,  2.2  and  2.3  were  all  developed  at  the  lattice  level  under  the 
assumption  of  isotropic,  homogeneous  material  properties  throughout  the  reference  volume  V.  This 
implies  that  the  the  weighted  exchange  energy  <l>o  employed  in  the  Helmholtz  relations  (4)  and 
(22)  is  constant  which  in  turn  implies  that  the  respective  bias  fields  Ej l  =  yjP-  and  H/ 1  =  yjj1 
and  parameters  a  =  and  a  =  respectively  employed  in  (17)  and  (25)  are  constant.  Since 
Ee  =  E+aP  and  He  =  H+aM  represent  effective  fields,  the  assumption  of  uniform  lattice  properties 
and  homogeneous  exchange  energy  relations  implies  constant  effective  fields  throughout  V. 

For  isotropic,  homogeneous,  single  crystal  materials,  the  local  average  models  (14),  (17)  or  (19) 
for  P,  (25),  (29)  or  (32)  for  M,  and  (43)  or  (50)  for  e  can  be  extended  throughout  the  compounds  to 
provide  bulk  or  macroscopic  constitutive  relations.  Hence  they  can  be  employed  to  quantify  single 
crystal  behavior  of  the  type  experimentally  measured  for  the  ferroelectric  material  BaTiOs  (e.g.,  see 
page  76  of  [72])  and  the  shape  memory  alloy  CuZnAl  [38,  109]. 

The  local  relations  do  not,  however,  accurately  quantify  the  gradual  transitions  and  pre-remanent 
switching  characteristic  of  nonhomogeneous,  polycrystalline  compounds  with  nonuniform  effective 
fields.  In  the  next  section,  we  employ  the  local  average  polarization,  magnetization  and  strain  rela¬ 
tions  as  kernels  when  developing  stochastic  homogenization  techniques  to  incorporate  the  variability 
inherent  to  the  majority  of  ferroic  compounds. 

3  Macroscopic  Models  for  Nonhomogeneous  Polycrystalline  Com¬ 
pounds 

The  mesoscopic  models  developed  in  Section  2  can  be  directly  employed  when  quantifying  the  polar¬ 
ization  and  stress-induced  strains  in  certain  single  crystal  ferroelectric  and  ferroelastic  compounds  as 
well  as  the  magnetization  in  some  uniaxial  wires  and  annealed  toroidal  specimens  [23] .  However,  the 
transitions  provided  by  these  local  models  are  too  steep  to  accurately  characterize  hysteresis  in  gen¬ 
eral  ferroic  materials  due  to  variations  in  the  free  energy  relations  produced  by  material,  stress  and 
field  nonhomogeneities,  nonunifornr  lattice  variations  across  grain  boundaries,  and  stress  and  crys¬ 
talline  anisotropies.  The  introduction  of  these  material  attributes  into  the  energy  relations  requires 
analysis  similar  to  that  employed  in  micromagnetic  models  [12,  13,  14,  32,  54]  or  micromechanical 
theory  [43]  and  typically  produces  models  whose  complexity  precludes  transducer  design  or  real¬ 
time  control  implementation.  Alternatively,  one  can  employ  local  energy  relations  as  kernels  from 
which  low-order  macroscopic  models  be  derived  either  through  homogenization  techniques  or  the 
determination  of  bulk  effective  parameters  through  stochastic  or  empirical  means.  In  this  section, 
we  consider  certain  parameters  in  the  local  energy-based  relations  to  be  stochastically  distributed 
to  reflect  variations  in  the  lattice,  grain  orientations,  or  exchange  energies.  This  yields  low-order 
macroscopic  models  in  which  the  majority  of  parameters  can  be  correlated  with  measured  properties 
of  experimental  data.  As  will  be  illustrated  in  Section  4,  the  models  accurately  characterize  both 
major  and  biased  minor  loop  data  for  a  wide  range  of  ferroic  compounds. 

3.1  Ferroelectric  Materials 

We  consider  first  the  construction  of  macroscopic  models  for  ferroelectric  compounds.  As  illustrated 
in  Figure  7,  variations  in  the  lattice  due  to  material  nonhomogeneities,  impurities,  grain  boundaries, 
or  polycrystallinity,  are  manifested  as  variation  in  the  Helmholtz  and  Gibbs  energies  which  in  turn 
produce  a  distribution  of  fundamental  hysterons  or  hysteresis  kernels.  To  incorporate  this  variability 
in  a  manner  which  facilitates  the  construction  of  low-order  macroscopic  models,  we  consider  the  local 
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Figure  7.  (a)  Nonuniform  lattice  and  polycrystalline  structure  of  PZT.  (b)  Helmholtz  energies 
associated  with  lattice  structures  (i)  and  (ii).  (c)  Variations  in  hysteresis  kernels  due  to  differing 
energy  profiles. 


coercive  fields  Ec,  given  by  (6),  to  be  manifestations  of  an  underlying  distribution  rather  than  fixed 
values  as  assumed  for  single  crystals  with  uniform  lattices.  To  enforce  the  physical  constraint  Ec  >  0, 
we  assume  that  variations  in  Ec  can  be  modeled  by  a  lognormal  distribution  with  the  density 

f(Ec)  =  Cie-MEJEc)/2c}i '  (53) 

It  is  noted  in  [114]  that  if  the  constant  c  is  small  compared  with  Ec,  the  mean  and  variance  of  the 
distribution  have  the  approximate  values 

(Ec)  «  Ec  ,  a  «  2 Ec  c .  (54) 

The  second  extension  of  the  local  average  relations  (14),  (18)  or  (21)  is  the  consideration  of 
effective  fields  in  the  material.  As  indicated  in  [6,  70],  the  applied  field  E  in  ferroelectric  materials  is 
augmented  by  fields  produced  by  neighboring  dipoles  which  produces  nonhonrogeneous  effective  fields 
Ee  in  the  materials.  This  produces  deviations  about  the  applied  field  which  can  produce  switching 
in  advance  of  the  renranence  point.  To  incorporate  these  variations,  we  consider  the  effective  field  to 
be  normally  distributed  about  the  applied  field  with  the  density 

f(Ee)  =  c2e-^E-E^/b.  (55) 
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The  macroscopic  polarization  model  combines  the  densities  (53)  and  (55)  to  yield 

poo  poo  _ 

[P(E)](t)  =  C  /  \P{Ee  +  E,Ec,0]{t)e-E2jbe-^E^E^/2c^dEedEc  (56) 

J  0  J —oo 

where  P  is  specified  by  (14),  (18)  or  (21)  and  £  denotes  the  initial  distribution  of  dipoles.  Discussion 
concerning  the  incorporation  of  lattice  variations  and  variable  effective  fields  in  the  Ising  relation 
(16)  as  well  as  details  regarding  the  approximation  of  the  integrals  and  implementation  of  the  model 
are  provided  in  [114]. 

We  note  that  while  the  polarization  model  (56)  does  incorporate  relaxation  mechanisms,  it  does 
not  incorporate  dynamic  elastic  effects  and  hence  should  be  employed  in  low  frequency  regimes. 
Initial  extensions  to  the  free  energy  relations  and  the  resulting  constitutive  relations  required  for 
constructing  fully  dynamic  models  have  been  investigated  in  [111].  Finally,  we  note  that  the  choice 
of  the  lognormal  and  normal  densities  (53)  and  (55)  is  based  on  their  mathematical  properties  rather 
than  a  priori  physical  arguments.  The  estimation  of  general  densities,  which  reflect  the  measured 
physical  properties  of  a  given  material  are  reported  in  [107]. 

3.2  Ferromagnetic  Materials 

To  accommodate  material  nonhomogeneities,  polycrystallinity  and  nonuniform  effective  fields  in 
ferromagnetic  materials,  we  employ  lattice  and  field  distributions  analogous  to  their  ferroelectric 
counterparts.  Lattice  variations  can  be  modeled  by  assuming  that  local  magnetic  coercive  fields  Hc 
are  manifestations  of  an  underlying  distribution  which,  for  the  models  presented  here,  we  assume  to 
have  the  lognormal  density 

f(Hc)  =  Cie-[HHc/Hc)/2c]\  (57) 

It  is  noted  in  [112],  where  the  present  modeling  framework  is  used  to  establish  an  energy  basis  for 
Preisach  models,  that  lognormal  densities  play  a  fundamental  role  in  Preisach  theory  for  magnetic 
materials  due  to  the  natural  manner  through  which  they  enforce  positivity  in  parameter  values  [31]. 

To  provide  additional  physical  motivation  regarding  assumed  effective  field  variability  and  to 
relate  this  analysis  to  previous  models,  we  consider  first  the  local  magnetization  relation  (25),  derived 
from  statistical  mechanics  tenets,  with  effective  fields  He  =  H  +  aM.  As  noted  in  Section  2,  the 
mean  field  parameter  a  is  given  by 

=  Hh  =  Ar4>0 
“  Ms  V  M| 

where  $0  is  the  energy  required  to  reorient  a  single  moment  in  an  ordered  lattice.  Furthermore,  it  is 
illustrated  in  [105]  that 

$o  =  2£J  (59) 

where  J  is  the  classical  exchange  integral  and  £  denotes  the  number  of  neighbors  adjacent  to  a  site. 

It  is  implicitly  assumed  in  the  mesoscopic  theory  yielding  M  that  J .  and  hence  $0,  is  constant. 
Thus  a  in  the  effective  field  He  =  H  +  aM  is  constant,  as  assumed  in  the  macroscopic  hysteresis 
models  developed  in  [56,  58].  However,  for  nonhomogeneous,  polycrystalline  materials,  it  is  more 
reasonable  to  assume  that  the  exchange  integral  J  is  variable  rather  than  constant.  At  the  quantum 
level,  this  variability  can  be  incorporated  by  modeling  the  overlap  of  electron  wave  functions  whereas 
at  the  macroscopic  level,  we  can  incorporate  this  variability  by  employing  effective  fields 

He  =  H  +  aM  (60) 
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where  a  satisfies  an  underlying  distribution.  To  avoid  the  implicit  dependence  on  M  afforded  by 
(60),  one  can  alternatively  model  or  estimate  the  distribution  of  He,  and  it  is  this  approach  that  we 
employ  here. 

The  assumption  that  the  effective  field  is  normally  distributed  about  the  applied  field  H  yields 
the  density 

7m  =  c2e-(H~H^b  (61) 


and  macroscopic  model 

/*oo  /'OO  _ 

[M(H)](t)  =  C  /  [M{He  + H,Hc,0\(t)e-H2/be-[ln(Hc/Hc)/2c^dHedHc 

J  0  J — oo 

/♦oo  /*oo 

=  c  /  v(Hc,He)\M(He  +  H,Hc,Z)\(t)dHedHc 

J  0  J — oo 


(62) 


for  the  bulk  magnetization  M.  The  kernel,  or  fundamental  hysteron,  M  is  given  by  (27),  (29)  or 
(32)  and  £  quantifies  the  initial  moment  orientation.  As  with  the  ferroelectric  model,  the  underlying 
assumption  of  lognormal  and  normal  densities  for  Hc  and  He  can  be  avoided  through  the  estimation 
of  a  general  density  v  through  techniques  analogous  to  those  employed  in  [102,  117]  for  estimating 
densities  in  Preisach  models. 

Details  regarding  the  formulation  and  implementation  of  (62)  can  be  found  in  [106,  105].  We  note 
that  because  the  given  formulation  does  not  incorporate  eddy  current  losses,  it  should  be  restricted 
to  low  frequency  regimes  or  transducer  configurations  which  minimize  eddy  current  losses. 


3.3  Ferroelastic  Materials 

The  effects  of  poly  crystallinity,  material  nonhomogeneities,  and  lattice  variations  across  grain  bound¬ 
aries  in  ferroelastic  materials  can  be  incorporated  by  considering  stochastic  distributions  of  the  ef¬ 
fective  stress  ae  and  either  the  transformation  stress  a  a  or  relative  stress  S  =  a  a  —  &m  given  by 
(42).  Because  cja  >  ffju,  as  depicted  in  Figure  6,  variability  in  6  can  be  incorporated  through  the 
assumption  that  5  is  is  lognormally  distributed  with  the  density  function 

f{6)  =  Cle-[In(5/5)/2c]2  (63) 

as  detailed  in  [69].  Furthermore,  it  is  illustrated  in  [97],  as  well  as  the  example  of  Section  4.3,  that 
for  certain  ferroelastic  compounds,  the  Laplace  relation 

f(aA)  =  (64) 

can  be  employed  to  accurately  characterize  the  distribution  of  a  a-  Two  appropriate  choices  for  the 
underlying  effective  field  distribution  are  the  normal  density  relation 

/Oe)  =  c2e-(CT-^)>  (65) 

or  the  Laplace  relation 

f(Ve)  =  (66) 

For  the  densities  (63)  and  (65),  this  yields  the  macroscopic  relation 

/*oo  /*oo  _ 

[e(a,  T)](t)  =  C  /  [e(ae  +  a,  5 ,  Ome^^e-^5^2 daed5  (67) 

J  0  J — oo 
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quantifying  bulk  strains  due  to  input  stresses  and  evolving  temperatures.  The  ferroelastic  kernel  e  is 
given  by  (43)  or  (50)  and  the  multiplicative  constant  is  simply  C  =  c\  ■  oi-  The  choice  of  the  Laplace 
densities  (64)  and  (66)  yields  the  analogous  relation 

POO  POO 

[e(a,T)](t)  =  C  /  [e(ae  + a,5,^)](t)e~^/be~^A~aA^cdaedaA  (68) 

J  0  J —oo 


for  e. 

3.4  Unified  Macroscopic  Models 

In  Section  2.4,  we  demonstrated  that  by  employing  a  general  order  parameter  e  =  P,M  or  e  and 
conjugate  field  <~p  =  E,  H  or  a,  mesoscopic  models  characterizing  hysteresis  in  ferroelectric,  ferromag¬ 
netic  and  ferroelastic  materials  could  be  consolidated  to  provide  a  unified,  energy-based  framework 
for  quantifying  hysteresis  at  the  lattice  level  in  ferroic  compounds.  Here  we  employ  these  lattice-level 
relations  as  kernels  to  construct  a  unified  macroscopic  framework  for  modeling  hysteresis  and  con¬ 
stitutive  nonlinearities  in  bulk  ferroic  materials.  This  framework  can  subsequently  be  employed  for 
material  characterization,  hybrid  transducer  design,  and  unified  control  design  for  systems  employing 
ferroelectric,  ferromagnetic  or  ferroelastic  actuators  and  sensors. 

We  let  (pc  =  Ec,  Hc  or  5  denote  general  coercive  or  relative  fields  and  let  ipe  =  Ee,  He  or  ae  denote 
effective  fields.  We  also  let  e  =  P,  M  or  e  denote  the  local  average  polarization,  magnetization,  and 
strain  respectively  defined  in  (14)  or  (21),  (27)  or  (32),  and  (43)  or  (50),  with  general  relations  for  e 
summarized  in  Section  2.4.  Initial  dipole,  moment  and  strain  configurations  are  specified  by  £. 

The  hysteresis  and  constitutive  nonlinearities  inherent  to  ferroic  compounds  can  then  be  quanti¬ 
fied  by  the  general  relation 

POO  POO 

[e(ip)\(t)  =  C  /  [e(</?e  +  <£,  Vc,0}{t)e-^2e/be~[hi^c/^c^2c]2dipedVc  (69) 

J  0  J — oo 

where  C,  b,Tpc  and  c  are  material-dependent  parameters.  If  additional  generality  is  required,  one  can 
employ  the  model 

POO  POO 

[e((p)\(t)  =  c  /  is(cpc,‘Pe)[e{‘Pe  + (70) 
J  0  J —oo 

where  v  is  a  general  density  to  be  estimated  through  a  least  squares  fit  to  data  or  direct  solution  tech¬ 
niques  analogous  to  those  developed  for  Preisach  models  [102].  As  indicated  in  previous  discussion, 
the  constitutive  models  incorporate  relaxation  mechanisms  but  do  not  incorporate  general  dynamic 
effects.  Hence  (69)  and  (70)  should  be  employed  in  low  frequency  regimes  or  modes  of  operation  for 
which  dynamic  effects  are  minimal. 

4  Experimental  Validation 

To  illustrate  the  performance  of  the  unified  modeling  framework,  we  consider  the  characterization  of 
hysteresis  and  constitutive  nonlinearities  in  PZT5A,  Terfenol-D,  and  NiTi.  In  all  cases,  quasistatic 
drive  conditions  were  employed  to  minimize  rate  dependencies  in  the  PZT  and  NiTi  and  eddy  current 
losses  in  the  Terfenol-D.  A  fundamental  step  in  the  model  construction  entails  the  estimation  of  model 
parameters  and  in  each  of  the  examples,  we  discuss  the  manner  through  which  initial  parameter 
estimates  can  be  obtained  using  properties  of  the  data. 


25 


4.1  Experimental  Validation  for  PZT5A 

We  consider  first  the  characterization  of  the  hysteretic  E-P  relation  for  PZT5A.  Data  was  collected 
from  a  1.7  cm  x  0.635  cm  x  0.0381  cm  PZT  wafer  at  200  mHz  with  peak  voltages  ranging  from  600  V 
to  1600  V  -  corresponding  field  values  can  be  computed  using  the  relation  E  =  V/h  where  h  =  3.81  x 
1CT4  nr  denotes  the  thickness  of  the  wafer  -  as  shown  in  Figure  8.  The  low  drive  frequency  yielded 
approximately  isothermal  operating  conditions  so  it  was  not  necessary  to  incorporate  temperature 
changes  through  the  evolution  equation  (15). 

The  field-polarization  relation  is  characterized  using  the  model  (56)  with  the  kernel  P  specified 
by  (21)  since  the  focus  here  is  on  the  quantification  of  multiple  drive  levels  rather  than  measurement 
and  characterization  of  relaxation  mechanisms.  We  point  out  that  similar  results  have  been  obtained 
with  the  kernel  (14)  which  retains  the  activation  energy  mechanisms  required  to  resolve  relaxation 
effects.  Regardless  of  kernel,  implementation  of  (56)  requires  the  approximation  of  the  integrals. 
This  was  accomplished  using  composite  Gaussian  quadrature  over  a  truncated  range  of  integration 
chosen  to  accommodate  the  decay  exhibited  by  the  densities  /  and  /.  As  illustrated  in  [114],  where 
details  regarding  the  implementation  of  the  model  are  provided,  convergence  was  obtained  for  both 
integrals  using  a  4  point  Gaussian  rule  over  a  composite  grid  having  20  divisions. 

To  construct  the  model,  the  parameters  Pr,  r],  b ,  Ec,  c  and  C  must  be  estimated  using  attributes  of 
the  data.  For  the  piecewise  linear  kernel  P  given  by  (21),  Pr  and  C  both  have  the  effect  of  scaling  the 
polarization  for  a  given  field  input.  Hence  they  can  be  combined  and  initially  chosen  to  obtain  correct 
saturation  behavior.  The  slope  of  the  kernel  scales  through  the  stochastic  homogenization  process 
so  measurements  of  the  reciprocal  slope  provide  initial  estimates  for  ?/.  From  (54),  it  follows  that 
an  initial  estimate  for  Ec  is  provided  by  the  measured  coercive  field  whereas  small  values  of  c  are 
required  when  characterizing  data  with  a  steep  transition  through  coercivity  since  little  variability  is 
exhibited  about  the  mean  value  Ec.  Finally,  the  parameter  b  quantifies  the  variance  in  the  effective 
field  which  determines,  in  part,  the  degree  to  which  switching  occurs  before  remanence.  Materials 


800  V 


Figure  8.  PZT5A  data  ( - )  and  model  predictions  ( - )  provided  by  (56)  with  the  kernel  P 

specified  by  (21).  Abscissas:  electric  field  (MV/m),  ordinates:  polarization  (C/m2). 
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with  nearly  linear  E-P  relations  at  renranence  yield  small  values  of  b  whereas  large  values  are  required 
to  accommodate  significant  switching  before  renranence. 

The  coercive  field  for  the  1600  V  data  is  1.2  x  106  V/nr  and  the  slope  after  reversal  at  saturation 
is  3.6  x  10s ,  and  these  two  values  were  used  as  initial  estimates  for  Ec  and  rj.  A  least  squares  fit  to  the 
1600  V  data  was  used  to  obtain  the  final  parameter  values  Pr  =  0.04  C/m2,  Ec  =  0.866010  xlO6  V/nr, 
rj  =  9.5  x  108,  c  =  0.4272  V2/nr2,  b  =  1.9754  x  1011  V2/nr2,  C  =  7.9926  x  10-12  yielding  the  model 
fit  illustrated  in  the  final  plot  of  Figure  8.  The  model  with  these  same  values  was  then  used  to 
predict  the  E-P  relations  for  600  V,  800  V  and  1000  V  inputs.  The  accuracy  of  both  the  high 
drive  level  fit  and  intermediate  drive  level  predictions  and  the  physical  nature  of  parameters  attests 
to  the  advantages  provided  by  combining  an  energy-based  kernel  with  stochastic  homogenization 
techniques  to  accommodate  material  nonhonrogeneities.  Further  examples  illustrating  the  predictive 
capabilities  of  the  model  for  PZT5H  and  PZT4  compounds,  as  well  as  numerical  examples  illustrating 
the  enforced  closure  of  multiply  nested  biased  minor  loops,  can  be  found  in  [114]. 

4.2  Experimental  Validation  for  a  Terfenol-D  Transducer 

To  illustrate  the  ferromagnetic  component  of  the  framework  as  well  as  the  capability  of  the  model 
to  guarantee  closure  of  biased  minor  loops  in  quasistatic  operating  regimes,  we  consider  the  char¬ 
acterization  H-M  relation  for  a  prototypical  Terfenol-D  transducer.  As  illustrated  in  Figure  9  and 
detailed  in  [17,  30],  Terfenol-D  transducers  employed  for  both  research  development  and  commercial 
applications  are  comprised  of  a  Terfenol-D  rod,  a  surrounding  wound  wire  solenoid,  a  biasing  perma¬ 
nent  magnet,  and  prestress  mechanisms.  Data  collected  at  0.2  Hz  under  zero  prestress,  isothermal 
room  temperature  conditions  is  plotted  in  Figure  10. 

The  hysteretic  and  nonlinear  relation  between  H  and  M  is  quantified  by  the  model  (62)  with 
the  kernel  M  specified  by  (27).  The  relationships  between  parameters  in  the  model  and  attributes 
of  the  data  are  analogous  to  the  ferroelectric  case  for  which  details  are  provided  in  Section  4.1. 
To  construct  the  model,  the  measured  coercive  field  value  Hc  =  6158  A/m  and  reciprocal  slope 

=  6.5  after  saturation  were  used  to  obtain  initial  estimates  for  Hc  and  p.  The  final  parameter 
values  Mji  =  8.7  x  104  A/m,  rj  =  7,  Hc  =  2000  A/m,  c  =  0.65  A2/m2,  b  =  5  x  108  A2/m2, 
C  =  1.98  x  10-8  were  then  obtained  through  a  least  squares  fit  to  the  symmetric  major  loop  data 
yielding  the  model  fit  illustrated  in  Figure  10.  Biased,  periodic  fields  were  subsequently  input  to 
the  model,  using  the  same  parameter  values,  to  obtain  the  minor  loop  predictions,  also  plotted  in 
Figure  10.  We  note  that  when  predicting  the  minor  loop  responses,  the  starting  magnetization  is 
determined  from  the  symmetric  loop  fit;  hence  the  accuracy  of  the  minor  loops  is  highly  dependent 
on  the  accuracy  of  the  symmetric  major  loop.  Additional  examples  illustrating  the  attributes  of  the 
hysteresis  model  for  ferromagnetic  materials  can  be  found  in  [105,  106]. 
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Figure  9.  Prototypical  Terfenol-D  transducer. 
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Figure  10.  Terfenol-D  data  and  model  response  provided  by  (62). 

4.3  Experimental  Validation  for  NiTi 

The  final  example  is  taken  from  [97]  and  illustrates  the  characterization  of  both  major  and  biased 
minor  loops  for  a  pseudoelastic  SMA  wire  in  a  tensile  experiment  conducted  at  room  temperature 
as  illustrated  in  Figure  11. 

The  wire  is  a  typical  NiTi  compound  (SE508)  from  Nitinol  Devices  and  Components  (NDC) 
having  a  diameter  of  d  =  0.5  mm  and  length  of  £  =  100  mm  which  underwent  a  training  procedure 
of  100  previous  loading/unloading  cycles.  The  mesoscopic  model  parameters  include  the  density 
p  =  6400  kg  nr-3,  the  two  Young’s  moduli  Y\  =  32  GPa  and  Ym  =  20  GPa,  and  the  transformation 
strain  et  =  0.038.  The  transformation  stresses  were  determined  to  be  a a{Tl)  =  430  MPa  at 
Ti  =  323  K  and  cta(Tu)  =  660  MPa  at  Tjj  =  353  K.  The  width  6  of  the  hysteresis  curve,  specified 
by  (42),  was  taken  to  be  5  =  100  MPa.  Since  convection  mechanisms  dominated  conduction  in  the 
experiments,  the  heat  transfer  coefficient  was  taken  to  be  hc  =  200  W  m~2  KG1  and  the  thermal 


Figure  11.  Superelastic  NiTi  data  [94]  and  model  fit  provided  by  (68)  with  e  specified  by  (50). 
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conductivity  was  specified  as  A  =  0  W  nr-1  K_1.  The  specific  heats  for  austenite  and  martensite 
were  taken  to  be  equal  and  having  the  value  of  C4  =  cm  =  450  J  kg^1  K_1.  The  relaxation  time  was 
chosen  to  be  T  =  10  ms  and  the  activation  volume  was  specified  to  be  1  x  10-23  nr3.  Finally,  the 
macroscopic  homogenization  was  performed  using  the  Laplace  distributions  (64)  and  (66)  with  the 
parameters  a  a  =  273  MPa  and  b  =  c  =  30  MPa. 

It  is  observed  that  the  model  accurately  characterizes  the  linearly  elastic  austenite  behavior,  the 
phase  transition  from  austenite  to  martensite  and  back  again,  and  the  behavior  of  the  nested  minor 
loops.  This  illustrates  the  utility  of  the  model  for  bulk  material  characterization,  transducer  design 
and  model-based  control  design. 

5  Concluding  Remarks 

The  theory  presented  here  provides  a  unified  framework  for  modeling  hysteresis  and  constitutive 
nonlinearities  inherent  to  a  broad  range  of  ferroelectric,  ferromagnetic,  and  ferroelastic  compounds. 
At  the  microscopic  scale,  the  sources  of  hysteresis  vary  quite  significantly  and  can  roughly  be  at¬ 
tributed  to  dipole  switching  in  ferroelectric  materials,  moment  rotation  and  domain  wall  losses  in 
ferromagnetic  materials,  and  transitions  between  austenite  and  martensite  phases  in  ferroelastic 
compounds.  At  the  lattice  level,  however,  analogous  Helmholtz  and  Gibbs  energy  relations  can  be 
formulated  and  employed  to  construct  unified  mesoscopic  models  both  in  the  presence  and  absence 
of  thermally  activated  relaxation  mechanisms.  For  homogeneous,  single  crystal  compounds,  the 
mesoscopic  relations  can  be  extrapolated  to  yield  macroscopic  hysteresis  models.  To  accommodate 
material  nonhomogeneities,  poly  crystallinity  and  nonunifornr  effective  fields,  stochastic  homogeniza¬ 
tion  techniques,  based  on  the  assumption  that  parameters  such  as  coercive,  relative  and  effective 
fields  are  manifestations  of  underlying  distributions,  are  employed  to  construct  macroscopic  models. 
These  models  provide  sufficient  accuracy  for  material  characterization  but  are  sufficiently  low-order 
to  permit  transducer  design  and  model-based  control  development. 

Given  the  general  nature  of  the  framework,  it  is  not  surprising  that  it  bears  certain  similarities 
with  Preisach  models  and  it  is  illustrated  in  [112]  that  it  actually  provides  an  energy  basis  for  ex¬ 
tended  Preisach  models  -  with  four  fundamental  differences,  (i)  Due  to  is  energy  basis,  a  number 
of  the  parameters  in  the  present  framework  can  be  directly  correlated  with  properties  of  the  data 
whereas  Preisach  parameters  must  typically  be  estimated  solely  through  a  least  squares  fit  to  data. 
While  least  squares  techniques  are  often  used  to  fine  tune  parameters  in  the  present  models,  the  phys¬ 
ical  correlation  provides  initial  estimates  and  facilitates  model  updating  to  accommodate  changing 
environmental  conditions,  (ii)  The  relaxation  mechanisms  and  temperature-dependence  are  incor¬ 
porated  in  the  kernels,  or  hysterons,  P,  M  and  e  rather  than  in  the  densities  or  parameters  as  is  the 
case  for  Preisach  models.  This  eliminates  the  necessity  of  vector-valued  parameters  or  lookup  tables 
in  subsequent  control  design,  (iii)  The  model  guarantees  biased  minor  loop  closure  in  quasistatic 
operating  regimes  but  does  not  enforce  congruency.  Hence  it  does  not  require  modifications  of  the 
type  employed  in  moving  Preisach  models  [31]  to  accommodate  the  noncongruency  often  encoun¬ 
tered  in  measured  data,  (iv)  The  model  automatically  incorporates  reversibility  through  the  energy 
construction  of  hysterons.  This  eliminates  the  modifications  required  in  Preisach  models  to  achieve 
reversibility  in  certain  drive  regimes. 

In  order  for  any  fixed-temperature,  rate  independent  model  to  provide  an  energy  basis  for  clas¬ 
sical  Preisach  representations,  it  must  satisfy  the  congruency  and  deletion  properties  established 
by  Mayergoyz  as  necessary  and  sufficient  conditions  for  classical  Preisach  models  [71].  Because  the 
energy-based  models  do  not  enforce  congruency,  they  are  thus  not  mathematically  equivalent  to 
classical  Preisach  formulations.  However,  they  do  provide  an  energy  basis  for  the  moving  Preisach 
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models  described  in  [31].  Moreover,  use  of  the  kernels  (14),  (32)  or  (50)  that  incorporate  relaxation 
mechanisms  through  statistical  mechanics  principles  yields  macroscopic  models  which  are  analogous 
to  the  statistically  modified  Preisach  models  used  to  model  aftereffects  and  accommodation  in  mag¬ 
netic  materials  [31].  Analysis  quantifying  the  relationship  between  this  modeling  framework  and 
various  extended  Preisach  theories  is  under  present  investigation. 

From  the  perspective  of  model-based  control  design,  this  modeling  approach  is  advantageous  due 
to  both  the  accuracy  and  efficiency  provided  by  the  models  and  the  fact  that  it  provides  a  unified 
framework  for  model-based  control  design.  In  the  context  of  nonlinear  optimal  control  design,  the 
models  have  recently  been  employed  when  developing  control  algorithms  for  piezoceramic  transduc¬ 
ers  [124]  and  SMA  actuators  [45,  95]  including  real-time  control  implementation  as  reported  in  [96]. 
To  permit  linear  control  design,  one  often  employs  inverse  models  as  filters  to  compensate  for  hys¬ 
teresis  and  constitutive  nonlinearities  so  that  control  inputs  to  the  system  are  approximately  those 
specified  by  the  control  law.  Whereas  a  number  of  control  designs  have  been  developed  within  this 
context  (e.g.,  see  [24,  26,  42,  44,  73,  74,  75,  103,  117,  118]),  all  are  dependent  on  having  accurate 
and  efficient  inverse  models.  For  homogeneous,  single  crystal  transducer  compounds,  the  kernels 
P ,  M  or  e  can  either  be  directly  inverted,  or  approximately  inverted  through  the  solution  of  comple¬ 
mentary  differential  equations  in  a  manner  similar  to  that  detailed  for  related  models  in  [103,  104]. 
For  nonhomogeneous,  polycrystalline  transducer  compounds,  approximate  model  inverses  have  been 
constructed  by  employing  the  highly  efficient  macroscopic  models  in  direct  algorithms  based  on  the 
monotonic  relation  between  external  fields  and  the  resulting  changes  in  P,  M  or  e  [111].  Robust 
control  designs  utilizing  approximate  inverses  constructed  in  this  manner  are  illustrated  for  rnagne- 
tostrictive  transducers  in  [74,  75],  and  analogous  inverse-based  control  designs  for  PZT  and  SMA 
transducers,  and  the  real-time  implementation  of  resulting  control  algorithms  are  under  present 
investigation. 
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